Numerical Simulation for Hydrogen-Assisted Cracking: An Explicit Phase-Field Formulation

Hydrogen-assisted cracking is one of the most dominant failure modes in metal hydrogen-facing materials. Therefore, the hydrogen-assisted cracking mechanism has been a hot topic for a long time. To date, there is very little published research on numerical methods to describe hydrogen-assisted cracking. This paper presents a new method for the description of hydrogen embrittlement crack growth: an explicit phase-field formulation, which is based on the phase-field description of cracks, Fick’s mass diffusion law, and the relationship between hydrogen content and fracture surface energy. A novel computational framework is then developed using the self-developed FEM software DYNA-WD. We numerically calculate several typical conditions in the 3-D coordinates to validate the effectiveness of the proposed computational framework. Specifically, we discuss (i) the failure of a square plate in a hydrogenous environment, (ii) the CT specimen failed with the inner hydrogen, (iii) the plate/failed with the corrosives, and (iv) the failure of the disk test. Finally, the relationship between Mises stress, the concentration of hydrogen, the thickness of the disc, and the loading rate is investigated.


Introduction
Hydrogen atoms will be able to permeate into the inner part of the material under the conditions of hydrogen, welding, surface heat treatment, and acid medium. Then, hydrogen atoms will be induced by stress to diffuse and enrich at the crack tip or stress concentration points. Under the action of high hydrogen concentration and high stress, the cracks spread quickly, resulting in what is known as hydrogen embrittlement at the low load [1][2][3][4][5]. Hydrogen embrittlement is a common problem in high-end manufacturing, transportation, and new energy industries. Especially in the new energy field, it is a key issue that needs to be urgently addressed during the storage, transport, and use process of high-pressure hydrogen containers.
In particular, the HELP and HEDE theories have been widely used in the explanation of the cracking of metallic materials. The former is based on the coupling of H2 and structural stress field, and the second is the effect of hydrogen concentration on the bonding

Method
In this section, we mainly discuss the computational framework of explicit phase-field formulation. As described above, the numerical scheme is involved in the field of phase-Materials 2023, 16, 1708 3 of 22 field method, Fick's mass diffusion law, FEM, and Euler method of differential equations. Focusing on these theories, we derived the calculation equations as follows.

Classical Phase-Field Fracture Model and Governing Equation
In the phase-field fracture model, on the finite element node, we first introduce a phasefield variable φ to describe the damage of the material. In one-dimensional coordinates, the phase-field variable φ, which is related to coordinate x, reads [3] φ(x) = e − |x−a| l 0 (1) here, a is the tip coordinate of the crack, and l 0 is the length-scale parameter, which represents the spread width of the crack. When l 0 tends to 0, φ represents a crack with almost zero width. In Equation (1), φ ∈ [0,1]. Moreover, φ = 0 and φ = 1 denote the intact and completely fragmented states of the material, respectively. The relationship between the distribution of φ and x, a, and l 0 is shown in Figure 1. cases. In addition, typical mixed-type failure, i.e., in the hydrogen embrittlement disk test, is numerically studied. Finally, concluding remarks are given in Section 4.

Method
In this section, we mainly discuss the computational framework of explicit phasefield formulation. As described above, the numerical scheme is involved in the field of phase-field method, Fick's mass diffusion law, FEM, and Euler method of differential equations. Focusing on these theories, we derived the calculation equations as follows.

Classical Phase-Field Fracture Model and Governing Equation
In the phase-field fracture model, on the finite element node, we first introduce a phase-field variable ϕ to describe the damage of the material. In one-dimensional coordinates, the phase-field variable ϕ, which is related to coordinate x, reads [3] (1) here, a is the tip coordinate of the crack, and l0 is the length-scale parameter, which represents the spread width of the crack. When l0 tends to 0, ϕ represents a crack with almost zero width. In Equation (1), ϕ ∈ [0,1]. Moreover, ϕ = 0 and ϕ = 1 denote the intact and completely fragmented states of the material, respectively. The relationship between the distribution of ϕ and x, a, and l0 is shown in Figure 1. Then, we compute the first and second derivatives ϕ′(x), ϕ″(x). The following relationship can be obtained analytically: Following [24][25][26][27], it is necessary to satisfy the boundary conditions in Equation (2), and the quadratic functional of the differential equation can be expressed as follows: Then, we compute the first and second derivatives φ (x), φ"(x). The following relationship can be obtained analytically: Following [24][25][26][27], it is necessary to satisfy the boundary conditions in Equation (2), and the quadratic functional of the differential equation can be expressed as follows: In 3-D coordinates, φ is a function of x, y, and z. At this time, the crack density function can be defined as [3,24,25] The gradient ∇φ can be written as During the transformation of a component from a continuum to a non-continuum, the fracture energy generated can be approximately calculated as [3] where G c (θ) is the energy release rate, Г is the crack surface, and Ω is the geometry.
In the continuum damage mechanics model, parameters to characterize material damage are introduced. Moreover, a kind of stress degradation function is introduced in the phase-field method. We couple it with the material constitutive calculation to describe the effect of cracks on components. Referring to [3], the stress degradation function can be defined as where a small amount k is introduced to avoid calculation termination caused by unit distortion. When calculating, k = 1 × 10 −7 [3]. When it does not take into account material damage, we define the strain energy as Ψ 0 ε . After material degradation, the strain energy can be calculated as Using the coupling effect of the stress field and the phase field, the function of the total potential energy can be expressed as Considering the existence of the stress degradation function g(Ψ ε ), the Cauchy stress can be obtained [21] as ∂Ψ 0 ε ∂ε (10) which is calculated by taking the first derivative of the strain energy with respect to the strain tensor. Accordingly, the governing equation based on phase-field fracture reads,

Hydrogen Diffusion Equation
In HELP and HEDE theories, it is found that hydrogen will diffuse under the stress field and phase field and reduce the aggregation energy between metal atoms. Thus, the material undergoes brittle fracture without undergoing a plastic phase. Alvaro et al. [28] performed extensive calculations based on first principles. Subsequently, they fitted the linear relationship between the critical energy release rate of the material and the hydrogen Materials 2023, 16, 1708 5 of 22 attachment rate θ. The result has been adopted by many scholars [1][2][3]21]; thus, the critical energy release rate G c (θ) at the hydrogen concentration θ can be defined as herein, G c (0) is the energy release rate of non-hydrogen-facing or non-hydrogen-containing materials, and χ is the hydrogen damage coefficient which indicates the degree of influence of hydrogen on the brittle fracture of materials. The surface hydrogen concentration θ can be calculated as where C is the volumetric hydrogen ion concentration, R is the universal gas constant, T is the temperature, and ∆g 0 b is Gibbs free energy. The diffusion of hydrogen in metal follows the law of conservation of mass, that is, it satisfies where ∂Ω q is the hydrogen concentration boundary. The above calculation must satisfy the boundary conditions. J is obtained as where D is the concentration diffusion coefficient, V H is the partial molar volume of hydrogen in the solid solution, and σ H is the hydrostatic pressure. Accordingly, the governing equation for hydrogen diffusion reads dC dt + ∇·J = 0 J·n = q on ∂Ω q (16) where q is the flux of hydrogen concentration. This equation satisfies the law of conservation of mass and the boundary conditions. In particular, it is noted that the Formulas (12)- (16) come from literature [3,23].

Spectral Decomposition of Strain Energy
After considering material damage, Equation (8) gives the calculation of strain energy. It is clear from this equation that the stress degradation function has an influence on the whole stress tensor. That is, in both tension and compression, energy is released when failure occurs. Of course, it is not physically realistic. Therefore, when the strain energy is not decomposed under tension, the analytical frame can only be used for the failure of the tensile stress, that is, the mode I failure. To overcome this problem, Miehe et al. [20] proposed a spectral decomposition method of strain energy. In 3-D coordinates, the tension and compression strain can be decomposed into [23] Herein, ε i is the main strain, <•> is Macaulay brackets. When damage is not considered, the strain energy can be calculated as where λ and µ are lame constants. After spectral decomposition, from the tensile and compressive strain energies generated by the tensile and compressive strains, we calculate, respectively, as Calculated by Equation (19), the strain energy after material degradation can be obtained as considering the material damage caused by tensile strain, we decompose the stress into tensile and compressive stress tensor. Tensile and compressive stress tensor can read as We note that Equations (18)-(21) are derived from [19].

C0 Shell Element
Referring to [29][30][31], the 3-node C0 shell element is shown in F where ̅( , ) is the position coordinate associated with natural co reference surface, the expression for which can be seen in Equation coordinate correction value along the thickness direction based on t vector [32].
where ̅ is the node coordinates and Na is the shape function at n The shape function at each node can be obtained as The Eulerian coordinates at its arbitrary point in the current configuration can be where x(ξ, η) is the position coordinate associated with natural coordinates (ξ, η) on the reference surface, the expression for which can be seen in Equation (23). X(ξ, η, ζ) is the coordinate correction value along the thickness direction based on the position coordinate vector [32].
where x a is the node coordinates and N a is the shape function at node a. The shape function at each node can be obtained as whereê a is the vector along the thickness direction at node a, and z a is the thickness direction function. The forces and torques at the nodes can be calculated as where B T M is the membrane strain and B T S is the contribution of the bending strain to the B matrix.

Finite Element Implementation of Computational Framework
Referring to [33], the displacement field u, the phase-field parameter φ, and the hydrogen concentration C at its arbitrary point in the cell can be calculated as, respectively, the following: In Equation (30), u is the generalized displacement. In the shell element, u = u x u y u zθxθy T . N u , N φ , and N C are the finite element shape functions [33].
The driving equations for the displacement field u, the phase-field parameter φ, and the hydrogen concentration C can be obtained as, in the explicit finite element calculation .
where u = {u e }, φ = {φ e }, and C = {C e }. Referring to [33], one can readily identify a mass matrix, a node external force, a node internal force, and C φ tensor, Accordingly, we obtain where (A) N e e=1 represents the element-to-global assembly in FEM, N e is the total number of units, N e u and N e φ both represent [N1, N2, N3], B e φ is the derivative of shape functions with respect to global coordinates [33,34]. ∇ · φ is the phase-field gradient, and ∇ · φ = B e φ φ.

Central Difference Method and Rodriguez Transformation
In the FEM explicit calculation, the plane displacement vector u p = [u x , u y , u z ], phasefield variable φ, and hydrogen concentration C can be updated by the central difference method [35][36][37]. The computational time axis is given in Figure 3. Here, t = 0 is the initial moment, t n−1/2 is the intermediate moment between t n−1/2 and t n , and ∆t n−1/2 is the difference between t n and t n−1 . Thus, the updated formula can be obtained as .
where ( ) represents the element-to-global assembly in FEM, is the total num of units, and both represent [ 1,2,3], is the derivative of shape fu tions with respect to global coordinates [33,34]. ∇ • is the phase-field gradient, and ∇ = .

Central Difference Method and Rodriguez Transformation
In the FEM explicit calculation, the plane displacement vector = [ , , ], pha field variable ϕ, and hydrogen concentration C can be updated by the central differe method [35][36][37]. The computational time axis is given in Figure 3. Here, t = 0 is the in moment, tn−1/2 is the intermediate moment between tn−1/2 and tn, and ∆tn−1/2 is the differe between tn and tn−1. Thus, the updated formula can be obtained as We use the Rodriguez transformation [38,39] to calculate rotation vector , shell elements, the angular coordinates at nodes are explicitly stored by the fiber vect Therefore, we use the Rodriguez transformation to update the fiber vector during ca lation. The vector rotation schematic is given in Figure 4. At t, we assume that the angu displacement vector of a certain node is e (t) = [e1 (t), e2 (t), e3 (t)], and the correspond angular velocity is w = [w1, w2, w3]. Accordingly, the vector e (t + ∆t) = [e1 (t + ∆t), e2 (t + e3 (t + ∆t)] at time t + ∆t can be calculated as The (∆ ) is given by We use the Rodriguez transformation [38,39] to calculate rotation vector θ x ,θ y . In shell elements, the angular coordinates at nodes are explicitly stored by the fiber vectors. Therefore, we use the Rodriguez transformation to update the fiber vector during calculation. The vector rotation schematic is given in Figure 4. At t, we assume that the angular displacement vector of a certain node is e (t) = [e 1 (t), e 2 (t), e 3 (t)], and the corresponding angular velocity is w = [w 1 , w 2 , w 3 ]. Accordingly, the vector e (t + ∆t) = [e 1 (t + ∆t), e 2 (t + ∆t), e 3 (t + ∆t)] at time t + ∆t can be calculated as

Numerical Experiment
From Equation (17), the most outstanding w tension and compression strain energy to the stre Therefore, it can be assumed that the original calcu I, mode II, and mixed failures. In this section, we an brittlement fracture for several typical thin-walled ment and internal hydrogen condition. Furthermo accuracy and effectiveness of our computational fra

Mode I Failure of Square Plate in a Hydrogenous E
As in [1][2][3]21], we study the crack propagatio hydrogenous environment and tensile loading, de plastic material model is used for the steel, and sing employed to discrete the square plate. Model geom are given in Figure 5. Firstly, we set the initial hy distributed and satisfy the following: The R ij (∆θ) is given by For ∆S ij , D m , and δ ij , where

Numerical Experiment
From Equation (17), the most outstanding work is considering the contribution of tension and compression strain energy to the stress field and phase field respectively. Therefore, it can be assumed that the original calculation framework can be used in mode I, mode II, and mixed failures. In this section, we analyze the calculation of hydrogen embrittlement fracture for several typical thin-walled structures in the hydrogen environment and internal hydrogen condition. Furthermore, we investigate the computational accuracy and effectiveness of our computational framework.

Mode I Failure of Square Plate in a Hydrogenous Environment
As in [1][2][3]21], we study the crack propagation process of a square steel plate in a hydrogenous environment and tensile loading, denoted as Case 1 for short. An elasticplastic material model is used for the steel, and single point integral C 0 shell elements are employed to discrete the square plate. Model geometry, loads, and boundary conditions are given in Figure 5. Firstly, we set the initial hydrogen concentration to be uniformly distributed and satisfy the following: C(t = 0) = C 0 [3]. Secondly, in the numerical calculation process, we set the boundary to be a constant hydrogen concentration, C = C b . All outer boundaries of the specimen are closely related to the environment, including the crack face. Thirdly, the parameter values are shown in Table 1, following [3,19].  We calculate the mechanical effects of square steel plates with different hydrogenfacing environments, that is, the values of C0 are 0, 0.1 wt ppm, 0.5 wt ppm, and 1.0 wt ppm. The calculated load-displacement curves are consistent with the results in reference [3]. Accordingly, the validity of the new computational framework has been validated. Figure 6 illustrates the load-displacement curves for different concentrations of hydrogen. It is anticipated that the greater the hydrogen concentration, the more vulnerable the member is to the hydrogen embrittlement, and the lower the critical failure load. It is worth noting that a finer grid size was used in our study. Compared with the coarse mesh, the critical failure load is slightly lower and the whole crack propagation lasts longer.
In this study, a displacement load in the tensile direction is applied on the upper boundary of the steel plate. As a result, the entire process is a model I failure. The stress is concentrated at the initial crack tip at the beginning. With the increase in the load, the crack starts from the tip and spreads horizontally. Based on the initial concentration of C0 = 0.5 wt ppm, the crack propagation process, i.e., the phase-field variable cloud diagram, is illustrated in Figure 7. The results obtained in this paper agree with those of [2,3,23].    In [3,23], the values of l 0 are 0.05 mm and 0.0075 mm, respectively. Their findings have shown that cracks show weaker accumulation due to the larger characteristic length. Correspondingly, the value is set to 0.002 mm in our study. To balance calculation accuracy and calculation efficiency, we use refined mesh in the failure area and stress concentration point, and coarse mesh in the rest. As shown in Figure 5b, the sizes of the refined and coarse meshes are 0.0025 and 0.015, respectively.
We calculate the mechanical effects of square steel plates with different hydrogenfacing environments, that is, the values of C 0 are 0, 0.1 wt ppm, 0.5 wt ppm, and 1.0 wt ppm. The calculated load-displacement curves are consistent with the results in reference [3]. Accordingly, the validity of the new computational framework has been validated. Figure 6 illustrates the load-displacement curves for different concentrations of hydrogen. It is anticipated that the greater the hydrogen concentration, the more vulnerable the member is to the hydrogen embrittlement, and the lower the critical failure load. It is worth noting that a finer grid size was used in our study. Compared with the coarse mesh, the critical failure load is slightly lower and the whole crack propagation lasts longer.
In this study, a displacement load in the tensile direction is applied on the upper boundary of the steel plate. As a result, the entire process is a model I failure. The stress is concentrated at the initial crack tip at the beginning. With the increase in the load, the crack starts from the tip and spreads horizontally. Based on the initial concentration of C 0 = 0.5 wt ppm, the crack propagation process, i.e., the phase-field variable cloud diagram, is illustrated in Figure 7. The results obtained in this paper agree with those of [2,3,23].
The length scale parameter l 0 determines the crack spread width. The distribution of the field variable φ after complete failure is shown in Figure 8, when l 0 takes values of 0.05 mm, 0.0075 mm, and 0.002 mm, respectively. As evident from Figure 8a, when the characteristic length is 0.05 mm, the calculated crack length is quite different from the actual one, and it is difficult to obtain accurate simulation results. When C 0 = 0.5 wt ppm, the load-displacement curves corresponding to different l 0 are shown in Figure 9. Obviously, when the length dimension is at the order of 10 −3 , the simulation results are in good agreement. However, the calculation results have a large deviation when the length dimension is at the order of 10 −2 (0.05 mm). Consequently, the length dimension is unsuitable for actual engineering.
worth noting that a finer grid size was used in our study. Compared with the the critical failure load is slightly lower and the whole crack propagation las In this study, a displacement load in the tensile direction is applied boundary of the steel plate. As a result, the entire process is a model I failu is concentrated at the initial crack tip at the beginning. With the increase in crack starts from the tip and spreads horizontally. Based on the initial conce = 0.5 wt ppm, the crack propagation process, i.e., the phase-field variable cl is illustrated in Figure 7. The results obtained in this paper agree with those Novel computing framework The length scale parameter l0 determines the crack spread width. The distribution of the field variable ϕ after complete failure is shown in Figure 8, when l0 takes values of 0.05 mm, 0.0075 mm, and 0.002 mm, respectively. As evident from Figure 8a, when the characteristic length is 0.05 mm, the calculated crack length is quite different from the actual one, and it is difficult to obtain accurate simulation results. When C0 = 0.5 wt ppm, the load-displacement curves corresponding to different l0 are shown in Figure 9. Obviously, when the length dimension is at the order of 10 −3 , the simulation results are in good agreement. However, the calculation results have a large deviation when the length dimension is at the order of 10 −2 (0.05 mm). Consequently, the length dimension is unsuitable for actual engineering. The length scale parameter l0 determines the crack spread width. The distribution of the field variable ϕ after complete failure is shown in Figure 8, when l0 takes values of 0.05 mm, 0.0075 mm, and 0.002 mm, respectively. As evident from Figure 8a, when the characteristic length is 0.05 mm, the calculated crack length is quite different from the actua one, and it is difficult to obtain accurate simulation results. When C0 = 0.5 wt ppm, the load-displacement curves corresponding to different l0 are shown in Figure 9. Obviously when the length dimension is at the order of 10 −3 , the simulation results are in good agreement. However, the calculation results have a large deviation when the length dimension is at the order of 10 −2 (0.05 mm). Consequently, the length dimension is unsuitable for actual engineering.     (14) and (15), the diffusion rate of hydrogen is positively correlated with the hydrostatic pressure σ H . Another significant finding is that stress concentration tends to occur at the crack tip. Therefore, during the stretching process, hydrogen ions diffuse to the crack tip [2,3]. As can be seen in Figure 10, hydrogen accumulates near the crack tip with increasing external load and hydrostatic stress. As is found in Equations (14) and (15), the diffusion rate of hydrogen is positiv correlated with the hydrostatic pressure σH. Another significant finding is that stress co centration tends to occur at the crack tip. Therefore, during the stretching process, hyd gen ions diffuse to the crack tip [2,3]. As can be seen in Figure 10, hydrogen accumula near the crack tip with increasing external load and hydrostatic stress.

As is found in Equations
With the increase in hydrogen concentration, the critical energy release rate is d creased. Therefore, this facilitates the initiation and subsequent propagation of cracks. F ure 6 illustrates the load-displacement curves at various concentrations of hydrogen. summary, our initial conclusions were confirmed.

Mode I Failure of CT Specimen with Internal Hydrogen
Residual hydrogen atoms are usually found in the local regions of the metal pa after casting, welding, electro-chemical processing, or surface heat treatment. In the ginning, the hydrogen atoms are caused by stress or stress concentration at the edge the crack, and then become rich. Furthermore, under the action of high hydrogen conce tration and high-stress interaction, the initiation and propagation of cracks occurred. last, hydrogen-assisted cracking takes place in the metallic elasticity phase. The pheno enon is known as internal hydrogen-assisted cracking. In response to this phenomeno Li et al. [40] conducted an experimental study on the cracking threshold KTH of hydrog embrittlement of high-performance martensitic steel AerMet100. Subsequently, Martín et al. carried out simulation studies corresponding to the experiments. Referring to [3] a [40], we investigate KTH, crack growth process, and load-displacement curves for A Met100 CT specimens, denoted as Case 2 for short. The geometry, loading, and finite e ment mesh of Case 2 are shown in Figure 11. Firstly, we employ the C0 shell to discr specimens, and we simulate the quasi-static loading with a 0.1 mm/min velocity load. S With the increase in hydrogen concentration, the critical energy release rate is decreased. Therefore, this facilitates the initiation and subsequent propagation of cracks. Figure 6 illustrates the load-displacement curves at various concentrations of hydrogen. In summary, our initial conclusions were confirmed.

Mode I Failure of CT Specimen with Internal Hydrogen
Residual hydrogen atoms are usually found in the local regions of the metal parts after casting, welding, electro-chemical processing, or surface heat treatment. In the beginning, the hydrogen atoms are caused by stress or stress concentration at the edge of the crack, and then become rich. Furthermore, under the action of high hydrogen concentration and high-stress interaction, the initiation and propagation of cracks occurred. At last, hydrogen-assisted cracking takes place in the metallic elasticity phase. The phenomenon is known as internal hydrogen-assisted cracking. In response to this phenomenon, Li et al. [40] conducted an experimental study on the cracking threshold K TH of hydrogen embrittlement of high-performance martensitic steel AerMet100. Subsequently, Martínez et al. carried out simulation studies corresponding to the experiments. Referring to [3,40], we investigate K TH , crack growth process, and load-displacement curves for AerMet100 CT specimens, denoted as Case 2 for short. The geometry, loading, and finite element mesh of Case 2 are shown in Figure 11. Firstly, we employ the C 0 shell to discrete specimens, and we simulate the quasi-static loading with a 0.1 mm/min velocity load. Secondly, the number of elements is 36,965 and the element size is 0.03 mm. Thirdly, Table 2 Figure 12 gives the crack growth process of the CT specimen with th ppm. The calculated results are in good agreement with those in reference [ Figure 12. Crack growth process of CT specimen.
Moreover, we compared the results of explicit phase-field formulatio  Note that K TH is an important parameter in the description of fracture mechanics failure. The results show that the crack will be unstable when the K TH reaches the critical value. The dimension of K TH is MPa/m 0.5 . Figure 12 gives the crack growth process of the CT specimen with the C 0 = 1.0 wt ppm. The calculated results are in good agreement with those in reference [3].   Figure 12 gives the crack growth process of the CT specimen with the C0 = 1.0 wt ppm. The calculated results are in good agreement with those in reference [3]. Moreover, we compared the results of explicit phase-field formulation for the fracture strength threshold KTH with experimental results in [40]. As can be seen from the results shown in Figure 13, the two particular analysis results are as follows.
1. The simulation results are in good agreement with the experimental results; 2. In the range of 0~1 ppm, KTH decreases sharply with the increase in hydrogen content.
In this case, the internal hydrogen embrittlement is especially severe for this kind of high strength steel, and at 1~8 ppm, it is not sensitive to the amount of hydrogen. Moreover, we compared the results of explicit phase-field formulation for the fracture strength threshold K TH with experimental results in [40]. As can be seen from the results shown in Figure 13, the two particular analysis results are as follows.
σy and the breaking strength when failure. It is considered that hydrogen reasons for the fracture of the CT, and at the concentration of hydrogen abo contribution of hydrogen to failure is essentially constant.

Mode I Failure for a Plate with Corrosion Pits
As is well known, metal failure is often caused by a combination of h other corrosion. For instance, in aqueous or acidic solutions, working com rapidly under the combined action of corrosion and hydrogen [42]. A commo is field pipes often rupturing from corrosion pits under the action of hydr other interesting occurrence is that castings and weldments often suffer fr embrittlement fractures from machining defects. Previous research has foun sion pits and defect points have complexity in shape and uncertainty in locat ber. Therefore, it is particularly necessary for the established calculation fram ble of describing the propagation of many cracks and the splitting and j cracks. In this section, we investigate the propagation of the crack under the displacement load on a sheet with three initial defects, which is referred to a K TH (Mpa·m -0.5 ) Figure 13. K TH as a function of the hydrogen concentration C in AerMet100 [40].

1.
The simulation results are in good agreement with the experimental results; 2.
In the range of 0~1 ppm, K TH decreases sharply with the increase in hydrogen content. In this case, the internal hydrogen embrittlement is especially severe for this kind of high strength steel, and at 1~8 ppm, it is not sensitive to the amount of hydrogen.
As is found in [41], Li et al. conducted a series of experimental studies on 30CrMo steel, and the second conclusion is similar to their experimental results. Figure 14 is a schematic drawing of a finite element with cracks in 2-D coordinates. In this case, the K TH can be represented as follows: ls 2023, 16,1708 σy and the breaking strength when failure. It is considered that h reasons for the fracture of the CT, and at the concentration of hydr contribution of hydrogen to failure is essentially constant.

Mode I Failure for a Plate with Corrosion Pits
As is well known, metal failure is often caused by a combin other corrosion. For instance, in aqueous or acidic solutions, wo Equation (53) shows that there is a positive correlation between K TH and the stress σ y , which is vertical to fracture direction. Due to hydrogen, there is a great difference between σ y and the breaking strength when failure. It is considered that hydrogen is one of the reasons for the fracture of the CT, and at the concentration of hydrogen above 1 ppm, the contribution of hydrogen to failure is essentially constant.

Mode I Failure for a Plate with Corrosion Pits
As is well known, metal failure is often caused by a combination of hydrogen and other corrosion. For instance, in aqueous or acidic solutions, working components fail rapidly under the combined action of corrosion and hydrogen [42]. A common occurrence is field pipes often rupturing from corrosion pits under the action of hydrogen [3]. Another interesting occurrence is that castings and weldments often suffer from hydrogen embrittlement fractures from machining defects. Previous research has found that corrosion pits and defect points have complexity in shape and uncertainty in location and number. Therefore, it is particularly necessary for the established calculation frame to be capable of describing the propagation of many cracks and the splitting and joining of the cracks. In this section, we investigate the propagation of the crack under the action of the displacement load on a sheet with three initial defects, which is referred to as Case 3. The geometry, loading, and FEM meshes of the defect are shown in Figure 15. ratio υ = 0.3, density ρ = 7900 kg/m 3 . At room temperature, the diffusion coefficient of iron D = 1 × 10 −8 mm 2 /s. The initial hydrogen concentration of the entire specimen is C0 = 1 wt ppm, and the loading speed u = 0.0416 μm/s. We mention that, even though the elastoplastic constitutive model has been applied to describe the mechanical performance of the sheet, it has no significance in setting the yield strength because of hydrogen embrittlement. The crack growth process is compared with that of literature [3], as shown in Figure  16. As can be seen from Figure 16, the error of the calculation results of the two methods is small. Significantly, it is close to the analytical result predicted by fracture mechanics in the literature [43]. The hydrogen diffusion mechanism on the meta-microscale has long been the focus of researchers [4]. On this case, it is concluded that the novel computational framework coupling RVE can be used to mode the hydrogen diffusion in the welding process.   In reference [3], the element size is set to 0.6 mm, and the characteristic length is 3.6 mm. In their simulation, the crack width is larger. Therefore, the characteristic length herein has a value of 0.8 mm. Besides, referring to [3], the following physical parameters are assumed: Young's modulus E = 20 GPa, energy release rate Gc(0) = 30 kJ/ m 2 , Poisson ratio υ = 0.3, density ρ = 7900 kg/m 3 . At room temperature, the diffusion coefficient of iron D = 1 × 10 −8 mm 2 /s. The initial hydrogen concentration of the entire specimen is C 0 = 1 wt ppm, and the loading speed u = 0.0416 µm/s. We mention that, even though the elastoplastic constitutive model has been applied to describe the mechanical performance of the sheet, it has no significance in setting the yield strength because of hydrogen embrittlement. The crack growth process is compared with that of literature [3], as shown in Figure 16. As can be seen from Figure 16, the error of the calculation results of the two methods is small. Significantly, it is close to the analytical result predicted by fracture mechanics in the literature [43].
The hydrogen diffusion mechanism on the meta-microscale has long been the focus of researchers [4]. On this case, it is concluded that the novel computational framework coupling RVE can be used to mode the hydrogen diffusion in the welding process.

Figure 16.
Comparison of crack growth process for Case 3 [3]. Figure 16. Comparison of crack growth process for Case 3 [3].

Mode II Failure for Disk Test
The disk test is mostly used to evaluate the hydrogen embrittlement susceptibility of metal materials. Compared with other tests, the disk test is more convenient to use, lower in cost, and higher in efficiency, due to simple equipment and small gas consumption. Accordingly, the disk test method is the most suitable method for the performance evaluation of hydrogen storage materials [44]. For standardization, the International Organization for Standardization (ISO) and the American Society for Testing Materials (ASTM) have formulated standards ISO 11114-4 and STM F1459-2006, respectively. Previous studies have demonstrated that the mechanical properties of metal materials in disk tests are related to strain rate [45], hydrogen concentration [46], and disk geometry. To verify the effectiveness of the proposed computational framework for mode II failure, we established the corresponding simulation model, denoted as Case 4 for short, referring to the tests in the literature [44]. The schematic diagram of the test equipment for Case 4 is shown in Figure 17.
According to the test, the geometric dimensions, constraint condition, and load of the disk in the model are given in Figure 18a. Following the literature [44], the test results show that the failure location is approximately within a circular area with a diameter of 20 mm. To improve the calculation efficiency, we make the mesh denser in the area, with a diameter of 30 mm. The mesh is given in Figure 18b. In the FE model, we set a total of 14,388 nodes and 28,690 elements. The minimum element size is 0.15 mm, and the characteristic length is 0.3 mm. We apply a full contraction on the blue region shown in Figure 18a. Additionally, the values of other physical parameters are as follows: Young's modulus E = 210 GPa, energy release rate Gc(0) = 30 kJ/m 2 , Poisson's ratio υ = 0.3, density ρ = 7900 kg/m 3 .
We numerically calculated the fracture process of high chromium alloy MANETA under the conditions of C 0 = 0, 1.5, 1.9, 2.5, 3.0, and 3.5 ppm. Moreover, we extracted the Mises stress at the crack tip when failure, denote by σ. The comparison between the calculated results and test results is shown in Table 3.
14388 nodes and 28690 elements. The minimum element size is 0.15 m teristic length is 0.3 mm. We apply a full contraction on the blue regio 18a. Additionally, the values of other physical parameters are as follow lus E = 210 GPa, energy release rate Gc(0) = 30 kJ/m 2 , Poisson's ratio υ 7900 kg/m 3 .   We numerically calculated the fracture process of high chromium alloy under the conditions of C0 = 0, 1.5, 1.9, 2.5, 3.0, and 3.5 ppm. Moreover, we ex Mises stress at the crack tip when failure, denote by . The comparison betwe culated results and test results is shown in Table 3.   Remark: Because of the manufacturing error of the disc, the thickness of the disc was not identical.
When C 0 = 3.5 ppm and the thickness of 0.525 mm, the propagation process of the hydrogen embrittlement crack is given in Figure 19. The hydrogen embrittlement cracks occur in a circular area with a radius of about 18 mm, which is consistent with the test results of 20 mm. It is important to note that the final failure pattern of the disc, as illustrated in Figure 20, is identical to that in the paper [44]. Remark: Because of the manufacturing error of the disc, the thickness of the disc was not identical When C0 = 3.5 ppm and the thickness of 0.525 mm, the propagation process of th hydrogen embrittlement crack is given in Figure 19. The hydrogen embrittlement crack occur in a circular area with a radius of about 18 mm, which is consistent with the tes results of 20 mm. It is important to note that the final failure pattern of the disc, as illus trated in Figure 20, is identical to that in the paper [44].  As shown in Table 3, the maximum error is 8.12% between numerical calculation and test results. So, the numerical computation frame can be applied to simulate the formation of hydrogen embrittlement cracks in 3-D for mode II failure.
When the diffusion speed reaches the speed of dislocation movement, hydrogen atoms segregate near the dislocation to hinder the movement of dislocation. Consequently, plastic loss of material is caused. It is reasonable to say that the strain rate is one of the most important control parameters in the disc test. Moreover, the thickness is one of the most important factors to influence the strength of thin-wall structures. Accordingly, we investigate the variation curve of , when C0 = 3.5 ppm, under different disk thicknesses h and different loading rates. The results are given in Figures 21 and 22 respectively.  Table 3, the maximum error is 8.12% between numerical calculation and test results. So, the numerical computation frame can be applied to simulate the formation of hydrogen embrittlement cracks in 3-D for mode II failure.

As shown in
When the diffusion speed reaches the speed of dislocation movement, hydrogen atoms segregate near the dislocation to hinder the movement of dislocation. Consequently, plastic loss of material is caused. It is reasonable to say that the strain rate is one of the most important control parameters in the disc test. Moreover, the thickness is one of the most important factors to influence the strength of thin-wall structures. Accordingly, we investigate the variation curve of σ, when C 0 = 3.5 ppm, under different disk thicknesses h and different loading rates. The results are given in Figures 21 and 22 respectively. When the diffusion speed reaches the speed of dislocation movement, hydrogen at-oms segregate near the dislocation to hinder the movement of dislocation. Consequently, plastic loss of material is caused. It is reasonable to say that the strain rate is one of the most important control parameters in the disc test. Moreover, the thickness is one of the most important factors to influence the strength of thin-wall structures. Accordingly, we investigate the variation curve of , when C0 = 3.5 ppm, under different disk thicknesses h and different loading rates. The results are given in Figures 21 and 22 respectively.  According to the 4th strength theory of material, the Mises press at one point can be written as follows: When failure occurs due to hydrogen, is less than σS, also known as the so-called hydrogen embrittlement. Smaller results in less hydrogen embrittlement resistance. In other words, the less there is, the more sensitive hydrogen embrittlement is. Figure 21 shows that resistance to hydrogen embrittlement is improved with the increase in the thickness h. In Figure 21, as the disk thickness h increases, also increases. It can be seen that increasing the thickness is still an effective means to improve the structural strength.
As outlined in Figure 22, when the loading rate is 10 MPa/min, the susceptibility to hydrogen embrittlement of the material is the highest. There is a more plausible explanation for this result. When the loading speed is low, such as 0.1 MPa/min, the loading process takes a long time. Consequently, the hydrogen atoms can fully diffuse to the crack tip so that the accumulated hydrogen atoms accelerate the plastic loss of the material. Comparatively, when the loading speed is large, such as 100 MPa/min, the large loading speed causes the whole loading process to not belong to quasi-static loading, so that the dynamic effect accelerates the crack growth. According to the 4th strength theory of material, the Mises press at one point σ can be written as follows: without hydrogen, when it meets Equation (55), this point is in yield phase.
When failure occurs due to hydrogen, σ is less than σ S , also known as the so-called hydrogen embrittlement. Smaller σ results in less hydrogen embrittlement resistance. In other words, the less σ there is, the more sensitive hydrogen embrittlement is. Figure 21 shows that resistance to hydrogen embrittlement is improved with the increase in the thickness h. In Figure 21, as the disk thickness h increases, σ also increases. It can be seen that increasing the thickness is still an effective means to improve the structural strength.
As outlined in Figure 22, when the loading rate is 10 MPa/min, the susceptibility to hydrogen embrittlement of the material is the highest. There is a more plausible explanation for this result. When the loading speed is low, such as 0.1 MPa/min, the loading process